In this script we are going to map Myeloid subtypes onto the Visium slides.
library(Seurat)
library(ggpubr)
library(cowplot)
library(dplyr)
library(ggplot2)
library(stringr)
library(readr)
library(SPOTlight)
library(SPATA2)
Loading necessary paths and parameters
source(here::here("misc/paths.R"))
source(here::here("utils/bin.R"))
"{myeloid}/{plt_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
"{myeloid}/{robj_dir}" %>%
glue::glue() %>%
here::here() %>%
dir.create(path = ,
showWarnings = FALSE,
recursive = TRUE)
set.seed(123)
Extract sample id and get Donor ID
# sample_id <- params$sample_id
sample_id <- "esvq52_nluss5"
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>% dplyr::pull(donor_id)
We have 8 different datasets that we are going to analyze separately. The spatial data comes from the script 03-clustering/03-clustering_integration.Rmd while the sc data can be found in Ramon’s scRNAseq analysis: /scratch/devel/rmassoni/tonsil_atlas_private/2-DOWNSTREAM_PROCESSING/results/R_objects/processed_seurat_objects/processed_seurat_objects/tonsil_integrated_with_harmony_scrublet_annotated.rds.
sp_obj <- "misc/{robj_dir}/20220215_tonsil_atlas_spatial_seurat_obj.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# Load SPOTlight data
spotlight_ls <- "{myeloid}/{robj_dir}/spotlight_ls_myeloid.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
# Load names df data
nm_df <- "{myeloid}/{robj_dir}/myeloid_nm_df.rds" %>%
glue::glue() %>%
here::here() %>%
readRDS(file = .)
Add colors to cell types
# (nm_df <- data.frame(col_vec))
# nm_df$plt_nm <- rownames(nm_df)
# nm_df$mod_nm <- stringr::str_replace_all(
# string = rownames(nm_df),
# pattern = "[[:punct:]]|[[:space:]]",
# replacement = ".")
decon_mtrx <- spotlight_ls[[2]]
decon_mtrx <- decon_mtrx[, colnames(decon_mtrx) != "res_ss"]
# Set as 0 cell types predicted to be under 3 % of the spot
decon_mtrx[decon_mtrx < 0.03] <- 0
Change column names
new_cn <- data.frame(mod_nm = colnames(decon_mtrx)) %>%
dplyr::left_join(nm_df, by = "mod_nm") %>%
dplyr::mutate(plt_nm = dplyr::if_else(is.na(plt_nm), mod_nm, plt_nm)) %>%
dplyr::distinct() %>%
dplyr::pull(plt_nm)
colnames(decon_mtrx) <- new_cn
We are going to add the deconvolution to the Seurat object.
sp_obj@meta.data <- cbind(sp_obj@meta.data, decon_mtrx)
Subset sample of interest
sp_sub <- sp_obj[, sp_obj@meta.data$gem_id == sample_id]
sp_sub@images <- sp_sub@images[Seurat::Images(sp_sub) == sample_id]
Check Topic profiles
nmf_mod_ls <- spotlight_ls[[1]]
nmf_mod <- nmf_mod_ls[[1]]
h <- NMF::coef(nmf_mod)
rownames(h) <- paste("Topic", 1:nrow(h), sep = "_")
topic_profile_plts <- SPOTlight::dot_plot_profiles_fun(
h = h,
train_cell_clust = nmf_mod_ls[[2]])
topic_profile_plts[[2]] +
ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90),
axis.text = ggplot2::element_text(size = 12))
Look at all cells profiles
topic_profile_plts[[1]] +
ggplot2::theme(
axis.text.y = ggplot2::element_blank(),
axis.text.x = ggplot2::element_blank(),
axis.title = ggplot2::element_blank())
Look at cells topic profile
basis_spotlight <- data.frame(NMF::basis(spotlight_ls[[1]][[1]]))
train_labs <- spotlight_ls[[1]][[2]]
colnames(basis_spotlight) <- unique(stringr::str_wrap(train_labs, width = 30))
basis_spotlight[basis_spotlight < 0.0000001] <- 0
DT::datatable(basis_spotlight, filter = "top")
Look at the location of each cell type in each slice separately
# Iterate over cell types
ct <- colnames(decon_mtrx)
# Iterate over images
lapply(names(sp_obj@images), function(nm) {
print(nm)
nm_donor <- id_sp_df %>% dplyr::filter(gem_id == nm) %>% dplyr::pull(donor_id)
# Iterate over cell types
ct_plt_ls <- lapply(ct, function(i) {
tmp_plt <- Seurat::SpatialFeaturePlot(
object = sp_obj,
features = i,
alpha = c(0, 1),
images = nm) +
ggplot2::scale_fill_gradientn(
colors = heat.colors(10, rev = TRUE)) +
ggplot2::scale_alpha(range = c(0, 1)) +
ggplot2::labs(title = stringr::str_wrap(string = i,
width = 25),
fill = "") +
ggplot2::theme(
plot.title = ggplot2::element_text(
hjust = 0.5,
size = 20,
face = "bold"))
if (sum(sp_sub@meta.data[, i]) == 0) {
tmp_plt <- suppressMessages(tmp_plt + ggplot2::scale_alpha(range = c(0,0)))
}
return(tmp_plt)
})
(plt_arr <- cowplot::plot_grid(
plotlist = ct_plt_ls,
axis = "trbl",
align = "hv",
nrow = 6,
ncol = 5))
"{myeloid}/{plt_dir}/cell_type_location_myeloid_{nm_donor}_new.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = plt_arr,
base_height = 30,
base_width = 25)
})
## [1] "tarwe1_xott6q"
## [1] "c28w2r_7jne4i"
## [1] "esvq52_nluss5"
## [1] "p7hv1g_tjgmyj"
## [1] "gcyl7c_cec61b"
## [1] "zrt7gl_lhyyar"
## [1] "qvwc8t_2vsr67"
## [1] "exvyh1_66caqq"
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-2-T_new.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-8-T_new.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-10-T_new.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-12-T_new.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-13-T_new.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-14-T_new.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-9-T_new.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cell_type_location_myeloid_BCLL-11-T_new.pdf"
Prepare data for boxplots
metadata_long <- sp_sub@meta.data %>%
# tidyr::pivot_longer(cols = c("annotation"),
# names_to = "stratification_id",
# values_to = "stratification_val") %>%
tidyr::pivot_longer(cols = all_of(ct), names_to = "ct_key", values_to = "ct_val") %>%
# dplyr::left_join(col_df, by = c("ct_key" = "ct_name")) %>%
dplyr::mutate(ct_val = dplyr::if_else(ct_val > 0.001, ct_val, 0))
Box plot of cell type proportion between stratified regions showing the unadjusted ANOVA Pvalue
keep_ct <- metadata_long %>%
dplyr::group_by(ct_key) %>%
dplyr::summarise(prop_sum = sum(ct_val)) %>%
dplyr::filter(prop_sum > 0) %>%
dplyr::pull(ct_key)
(bplt <- metadata_long %>%
dplyr::filter(ct_key %in% keep_ct) %>%
dplyr::mutate(
ct_key = stringr::str_wrap(string = ct_key,
width = 30)) %>%
ggpubr::ggboxplot(data = .,
x = "annotation_20220215",
y = "ct_val",
facet.by = "ct_key",
color = "annotation_20220215",
fill = "annotation_20220215",
add = "jitter",
scales = "free",
repel = TRUE,
outlier.shape = NA,
alpha = 0.6,
palette = "Set1",
ncol = 6) +
ggplot2::theme(
strip.text = ggplot2::element_text(size = 18, face = "bold"),
axis.text.y = ggplot2::element_text(size = 16),
axis.title.y = ggplot2::element_text(size = 20),
# axis.text.x = element_text(size = 12, angle = 90, vjust = 0.5, hjust = 0.5),
axis.text.x = ggplot2::element_blank(),
legend.text = ggplot2::element_text(size = 18),
legend.title = ggplot2::element_blank(),
strip.background = ggplot2::element_blank()) +
ggplot2::labs(
y = "Proportion",
color = "Regions",
fill = "Regions"))
# bplt <- bplt +
# ggpubr::stat_compare_means(method = "anova", size = 6) +
# ggplot2::scale_y_continuous(
# expand = expansion(mult = c(0, 0.1)),
# labels = function(x) sprintf("%.2f", x))
"{myeloid}/{plt_dir}/strat_bplot_{donor_id}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = bplt,
base_height = 20,
base_width = 25)
Slide-specific
lapply(Images(sp_obj), function(i) {
se_sub <- subset(sp_obj, subset = gem_id == i)
# se_sub
se_sub@images <- se_sub@images[Seurat::Images(se_sub) == i]
donor_id <- id_sp_df[id_sp_df$gem_id == sample_id, ] %>%
dplyr::pull(donor_id)
(cor_mtrx_ct <- SCrafty::correlation_heatmap(
se = sp_obj,
feats = ct,
assay = "Spatial",
slot = "data") +
ggplot2::labs(
title = glue::glue("{donor_id}: Cell-type correlation matrix")))
"{myeloid}/{plt_dir}/cor-mtrx_cell-type_{donor_id}.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_mtrx_ct,
base_height = 9,
base_width = 10)
})
## [[1]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[2]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[3]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[4]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[5]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[6]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[7]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
##
## [[8]]
## [1] "/scratch/devel/melosua/phd/projects/BCLLatlas/tonsil_atlas/spatial_transcriptomics/myeloid_integration/2020-09-22/plots_2020-09-22/cor-mtrx_cell-type_BCLL-10-T.pdf"
Integrated
# se_sub <- subset(merged_se, subset = gem_id == "esvq52_nluss5")
# se_sub
# se_sub@images <- se_sub@images[Seurat::Images(se_sub) == "esvq52_nluss5"]
(cor_mtrx_ct <- SCrafty::correlation_heatmap(
se = sp_obj,
feats = ct,
assay = "Spatial",
slot = "data") +
ggplot2::labs(
title = "Integrated cell-type correlation matrix"))
"{myeloid}/{plt_dir}/cor-mtrx_cell-type_integrated.pdf" %>%
glue::glue() %>%
here::here() %>%
cowplot::save_plot(
filename = .,
plot = cor_mtrx_ct,
base_height = 9,
base_width = 10)
sessionInfo()
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-conda_cos6-linux-gnu (64-bit)
## Running under: Red Hat Enterprise Linux Server release 6.7 (Santiago)
##
## Matrix products: default
## BLAS/LAPACK: /scratch/groups/singlecell/software/anaconda3/envs/spatial_r/lib/libopenblasp-r0.3.12.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C LC_TIME=es_ES.UTF-8 LC_COLLATE=en_US.UTF-8 LC_MONETARY=es_ES.UTF-8 LC_MESSAGES=en_US.UTF-8 LC_PAPER=es_ES.UTF-8 LC_NAME=C LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] parallel stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidyr_1.2.0 tibble_3.1.6 Biobase_2.50.0 BiocGenerics_0.36.1 SPATA2_0.1.0 SPOTlight_0.1.7 readr_2.1.2 stringr_1.4.0 dplyr_1.0.8 cowplot_1.1.1 ggpubr_0.4.0 ggplot2_3.3.5 SeuratObject_4.0.4 Seurat_4.1.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.2 reticulate_1.24 tidyselect_1.1.2 htmlwidgets_1.5.4 grid_4.0.1 Rtsne_0.15 munsell_0.5.0 codetools_0.2-18 ragg_1.2.2 ica_1.0-2 DT_0.21 future_1.24.0 miniUI_0.1.1.1 withr_2.5.0 spatstat.random_2.1-0 colorspace_2.0-3 highr_0.9 knitr_1.37 stats4_4.0.1 SingleCellExperiment_1.12.0 ROCR_1.0-11 ggsignif_0.6.3 tensor_1.5 listenv_0.8.0 NMF_0.23.0 MatrixGenerics_1.2.1 labeling_0.4.2 GenomeInfoDbData_1.2.4 polyclip_1.10-0 bit64_4.0.5 farver_2.1.0 rprojroot_2.0.2 parallelly_1.30.0 vctrs_0.3.8 generics_0.1.2 xfun_0.30 R6_2.5.1 doParallel_1.0.17 GenomeInfoDb_1.26.7 bitops_1.0-7 spatstat.utils_2.3-0 DelayedArray_0.16.3
## [43] assertthat_0.2.1 promises_1.2.0.1 scales_1.1.1 vroom_1.5.7 gtable_0.3.0 globals_0.14.0 goftest_1.2-3 rlang_1.0.2 systemfonts_1.0.4 splines_4.0.1 rstatix_0.7.0 lazyeval_0.2.2 spatstat.geom_2.3-2 broom_0.7.12 yaml_2.3.5 reshape2_1.4.4 abind_1.4-5 crosstalk_1.2.0 backports_1.4.1 httpuv_1.6.5 tools_4.0.1 gridBase_0.4-7 ellipsis_0.3.2 spatstat.core_2.4-0 jquerylib_0.1.4 RColorBrewer_1.1-2 ggridges_0.5.3 Rcpp_1.0.8 plyr_1.8.6 zlibbioc_1.36.0 purrr_0.3.4 RCurl_1.98-1.6 rpart_4.1.16 deldir_1.0-6 pbapply_1.5-0 S4Vectors_0.28.1 zoo_1.8-9 SummarizedExperiment_1.20.0 ggrepel_0.9.1 cluster_2.1.2 here_1.0.1 magrittr_2.0.2
## [85] data.table_1.14.2 scattermore_0.8 lmtest_0.9-39 RANN_2.6.1 fitdistrplus_1.1-6 matrixStats_0.61.0 hms_1.1.1 patchwork_1.1.1 mime_0.12 evaluate_0.15 xtable_1.8-4 IRanges_2.24.1 gridExtra_2.3 compiler_4.0.1 KernSmooth_2.23-20 crayon_1.5.0 htmltools_0.5.2 mgcv_1.8-39 later_1.3.0 tzdb_0.2.0 DBI_1.1.2 corrplot_0.92 MASS_7.3-55 Matrix_1.4-0 car_3.0-12 cli_3.2.0 SCrafty_0.1.0 igraph_1.2.11 GenomicRanges_1.42.0 pkgconfig_2.0.3 registry_0.5-1 plotly_4.10.0 spatstat.sparse_2.1-0 foreach_1.5.2 ggcorrplot_0.1.3 bslib_0.3.1 rngtools_1.5.2 pkgmaker_0.32.2 XVector_0.30.0 digest_0.6.29 sctransform_0.3.3 RcppAnnoy_0.0.19
## [127] spatstat.data_2.1-2 rmarkdown_2.12 leiden_0.3.9 uwot_0.1.11 shiny_1.7.1 lifecycle_1.0.1 nlme_3.1-155 jsonlite_1.8.0 carData_3.0-5 viridisLite_0.4.0 fansi_1.0.2 pillar_1.7.0 lattice_0.20-45 fastmap_1.1.0 httr_1.4.2 survival_3.3-1 glue_1.6.2 png_0.1-7 iterators_1.0.14 bit_4.0.4 stringi_1.7.6 sass_0.4.0 textshaping_0.3.6 irlba_2.3.5 future.apply_1.8.1